In [1]:
using QuadGK, Plots, Roots, ForwardDiff, Interpolations, LaTeXStrings, Polynomials, PlotlyJS
using SpecialFunctions, LinearAlgebra
xb=10
RR=5
Gamma = 1.22542
masa=100
rho0 = masa/(2 *sqrt(2)*1.22542*pi*RR^3 )*0.1 ### zmienione zostalo o 0.05
G=1
gamma=0.23
delta=4*sqrt(3)*pi*0.23


plotlyjs()
#gr()
function energydensityNH(x0)
return rho0*exp(-x0^2/(2*RR^2))
end
    
function energydensityOS(x0)

if x0<=10        
ret = rho0
        elseif x0>10
ret = 0
end
return ret
end            
            
#function Beta0(x0)
#return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * quadgk(x->x^2*rho0*exp(-x^2/(2*RR^2)),0,x0)[1] ) )
#end
function Beta0(x0)
    
if x0<=10
return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * (x0^3)/3*rho0 ) )

        
    elseif x0>10
    return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * (10^3)/3*rho0 ) )
        
    end
#return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * quadgk(x->x^2*energydensityOS(x),0,x0)[1] ) )
    
    #return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * quadgk(x->x^2*rho0*exp(-x^2/(2*RR^2)),0,x0)[1] ) )
#return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * rho0*(-exp(-x0^2/(2*RR^2))*RR^2*x0 +sqrt(pi/2)*RR^3*erf(x0/(sqrt(2)*RR ))   ) ) )
end

function atan2(xx)
if xx<0
return atan(xx) +pi
    else 
        return atan(xx) 
    end
end




function beta(lam,beta0)
    
return -atan2(1/(1/2*(-3*lam-2*1/tan(beta0)) ))
end


function x(lam, x0)
beta0=Beta0(x0)
return (x0*(4+9*lam^2 + 12*lam*1/tan(beta0) +4*1/(tan(beta0))^2  )^(1/3))/(2^(2/3) * (1+ 1/(tan(beta0))^2 )^(1/3) )
end

function t(lam)
return gamma*sqrt(delta)*lam
end

function lambda(t)
return t/(gamma*sqrt(delta))
end

function szift(lam, x0)
return -x(lam,x0)/(gamma*sqrt(delta))*sin(beta(lam,Beta0(x0)))*cos(beta(lam,Beta0(x0)))
end

#function xprim(lam, x0)
#beta0=Beta0(x0)
#function xx(lam)
#return (x0*(4+9*lam^2 + 12*lam*1/tan(beta0) +4*1/(tan(beta0))^2  )^(1/3))/(2^(2/3) * (1+ 1/(tan(beta0))^2 )^(1/3) )
#end
    
#xp = ForwardDiff.derivative(xx,lam)    
    
#return xp
#end



function xprim(lam, x0)
return x(lam, x0)*cos(beta(lam,Beta0(x0)))*sin(beta(lam,Beta0(x0))) 
end   

#Plots.plot(range(-30,30,length=1000),[lam->x(lam,15),lam->xprim(lam,15)])

function tprim(lam)
return gamma*sqrt(delta)
end


Plots.plot(range(-100,100,length=1000), lam -> beta(lam, Beta0(10)) )

The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.

Out[1]:
In [2]:
p=Plots.plot()
for x01 in range(0,30,length=30)
Plots.plot!(range(0,100,length=500),l->x(l,x01))
end
display(p)
In [3]:
x01=10
Plots.plot(range(0,100,length=500),l->x(l,x01), color=:red)

x01=5
Plots.plot!(range(0,100,length=500),l->x(l,x01), color=:blue)


x01=11
Plots.plot!(range(0,100,length=500),l->x(l,x01))

x01=15
Plots.plot!(range(0,100,length=500),l->x(l,x01))
Out[3]:
In [4]:
function fff(lam, x0)
return -1/(1-szift(lam,x0))
end

function derfpox0(lam,x0)

function test(x0)
return fff(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end



function derxpox0(lam, x0)
    
function test(x0)
return x(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end  


function derxpox0(lam, x0)

gamma1=0.23
rr=RR   
x0=x0+im*0  
lam=lam+im*0
rr=RR

ret=(-1 * rr ^ -1 * x0 ^ 2 * ((-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1) ^ (2//3) * ((2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (-1 * (x0 ^ 3 + 9 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * lam ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1)) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 6 * gamma1 * lam * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-(1//2) * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 18 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * lam ^ 2 * rr ^ 2)) ^ (-2//3) * (G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) ^ (-1//2) * (x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) ^ (-1//2) * (-1 * rr * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-(1//2) * rr ^ -2 * x0 ^ 2)) * (6 * pi * G * delta * rho0 * gamma1 ^ 2 * lam ^ 2 + exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp((1//4) * rr ^ -2 * x0 ^ 2) + 2 * G * delta * gamma1 * lam * rho0 * x0 * sqrt(pi) * ((x0 ^ 2 + -3 * rr ^ 2) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 16 * pi * G * delta * rho0 * gamma1 ^ 2 * rr ^ 2) + pi * G * delta * gamma1 * lam * rho0 * sqrt(2) * rr ^ 3 * (3 * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -16 * pi * G * delta * rho0 * gamma1 ^ 2) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp(-(5//12) * rr ^ -2 * x0 ^ 2))

return real(ret) 
end  


function derPRAWOpolam(lam, x0)
function test(lam)
return fff(lam,x0)*derxpox0(lam, x0)
end
return ForwardDiff.derivative(test, lam)
end  

#@manipulate for lam1 in range(-40,40,length=300)


#Plots.plot(range(0.1,40,length=1000), x0-> gamma*sqrt(delta)*derfpox0(lam1,x0)+derPRAWOpolam(lam1, x0), ylims=(-10,10))
#end
Out[4]:
derPRAWOpolam (generic function with 1 method)
In [5]:
function x0odx(x0, lam)
    


lam1=lam

tab1=range(0.001,100,length=1000)
tab2=x.(lam1,tab1)

xx=linear_interpolation(tab2,tab1)

return xx(x0)    
end


x0odx(11, 0)


#@manipulate for lam1 in range(-40,40,length=100)
##
#
#Plots.plot(range(0.1,40,length=1000), x1-> gamma*sqrt(delta)*derfpox0(lam1,x0odx(x1,lam1))
#        +derPRAWOpolam(lam1, x0odx(x1,lam1)), ylims=(-10,10))
#end
Out[5]:
10.999999999999998
In [6]:
function derxpox0polam(lam, x0)
function test(lam)
return derxpox0(lam, x0)
end
return ForwardDiff.derivative(test, lam)
end   


function derxpox0polam(lam, x0)

gamma1=0.23
rr=RR   
x0=x0+im*0  
lam=lam+im*0
rr=RR

ret=(-1 * rr ^ -1 * x0 ^ 2 * ((-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1) ^ (2//3) * ((2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (-1 * (x0 ^ 3 + 9 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * lam ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1)) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 6 * gamma1 * lam * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 18 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * lam ^ 2 * rr ^ 2)) ^ (-2//3) * (G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) ^( -1//2) * (x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^( 3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) ^ (-1//2) * (2 * G * delta * gamma1 * rho0 * x0 * sqrt(pi) * ((x0 ^ 2 + -3 * rr ^ 2) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 16 * pi * G * delta * rho0 * gamma1 ^ 2 * rr ^ 2) + pi * G * delta * gamma1 * rho0 * sqrt(2) * rr ^ 3 * (3 * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -16 * pi * G * delta * rho0 * gamma1 ^ 2) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -12 * pi * G * delta * lam * rho0 * rr * gamma1 ^ 2 * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^( 3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2)) * exp(-5//12 * rr ^ -2 * x0 ^ 2) + (2//3) * rr ^ -1 * x0 ^ 2 * ((-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1) ^( 2//3) * ((2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (-1 * (x0 ^ 3 + 9 * G * delta * rho0 * sqrt(2) * pi ^( 3//2) * gamma1 ^ 2 * lam ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1)) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 6 * gamma1 * lam * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 18 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * lam ^ 2 * rr ^ 2)) ^ (-5//3) * (G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) ^ (-1//2) * (2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) ^( -1//2 ) * (-1 * rr * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * (6 * pi * G * delta * rho0 * gamma1 ^ 2 * lam ^ 2 + exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 2 * G * delta * gamma1 * lam * rho0 * x0 * sqrt(pi) * ((x0 ^ 2 + -3 * rr ^ 2) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 16 * pi * G * delta * rho0 * gamma1 ^ 2 * rr ^ 2) + pi * G * delta * gamma1 * lam * rho0 * sqrt(2) * rr ^ 3 * (3 * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -16 * pi * G * delta * rho0 * gamma1 ^ 2) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * (6 * gamma1 * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 36 * pi * G * delta * lam * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 + -18 * G * delta * lam * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp(-5//12 * rr ^ -2 * x0 ^ 2))
return real(ret)    
end
Out[6]:
derxpox0polam (generic function with 1 method)
In [7]:
using DifferentialEquations

xpoczU=10

function rozwiazanie(s1,lamstar)

if s1<0    
    
function calyproblem(du, u, p, s)
du[1] = -derxpox0(u[1], u[2])
du[2] = -gamma*sqrt(delta)
du[3] = u[3]*derxpox0polam(u[1], u[2])
end

u0 = [lamstar; xpoczU; -1/(1-szift(lamstar,xpoczU))]

sspan=(0,-s1 ) 

prob3 = ODEProblem(calyproblem, u0, sspan, abstol = 1e-5, reltol = 1e-5)
sol3 = solve(prob3)

return sol3(-s1)
        
    elseif s1>=0        
        
function calyproblem2(du, u, p, s)
du[1] = derxpox0(u[1], u[2])
du[2] = gamma*sqrt(delta)
du[3] = -u[3]*derxpox0polam(u[1], u[2])
end

u02 = [lamstar; xpoczU; -1/(1-szift(lamstar,xpoczU))]

sspan2=(0,s1+1 ) 

prob32 = ODEProblem(calyproblem2, u02, sspan2, abstol = 1e-5, reltol = 1e-5)
sol32 = solve(prob32)

return sol32(s1)        
        
end
end

plotlyjs()
p=Plots.plot()
for i1 in range(-40,140,length=30)

    
Plots.plot!( s-> rozwiazanie(s,i1)[2], s->rozwiazanie(s,i1)[1], range(-20,10,length=500) )
end

display(p)
#Plots.plot(range(1e-15,10,length=1000),s-> sol3(s)[3])
In [8]:
function lambda1(s,lamstar)
return rozwiazanie(s,lamstar)[1]
end

function sodx0(x0)
return (x0-xpoczU)/(gamma*sqrt(delta))
end


function lamstar1(x0, lambda)
    
s=sodx0(x0)    
tab1=range(-150,150,length=1000)
tab2=lambda1.(s,tab1)

itp=linear_interpolation(tab2,tab1)

return itp(lambda)
end


function lamstar2(x0, lambda)
    
s=sodx0(x0)    

function test(l)    
return lambda1(s,l)-lambda
end     

ret=find_zero(test,-400)

return ret
end

#lamstar1(1,lambda1(-5,21))

x001=0.05
lamstar1(x001,lambda1(sodx0(x001),100))

function Ffinal(lam, x0)
    
if x0<xpoczU
    
ret = rozwiazanie(sodx0(x0),lamstar2(x0,lam))[3]
return ret
        
        
    elseif x0>=xpoczU       
return -1/(1-szift(lam,x0)) 
            
end 
        
        
end

#Plots.plot(range(-15,35,length=500), lamstar->lambda1(sodx0(1),lamstar) )
Out[8]:
Ffinal (generic function with 1 method)
In [9]:
function derFpox0(lam,x0)

function test(x0)
return Ffinal(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end


function derFpox02(lam,x0)

function test(x0)
return Ffinal(lam, x0)
end
    
dx0=1e-6
    
return (test(x0+dx0)-test(x0))/dx0
end



function derF_PRAWOpolam(lam, x0)
    
    
function test(lam)
return Ffinal(lam,x0)*derxpox0(lam, x0)
end
    
dlam=1e-6
    
return (test(lam+dlam)-test(lam))/dlam
end  


#@time derFpox02(1,5)

#@time derF_PRAWOpolam(1, 1)

#lam1=20
#@time Plots.plot(range(0.1,40,length=500), x0-> gamma*sqrt(delta)*derFpox02(lam1,x0)+derF_PRAWOpolam(lam1, x0), ylims=(-10,10))
Out[9]:
derF_PRAWOpolam (generic function with 1 method)
In [10]:
function lampos(s,lamstar)
    
function test(s)    
return rozwiazanie(s,lamstar)[1]
end
    
return ForwardDiff.derivative(test,s)
end

function lampos2(s,lamstar)

ds=1e-7
ret=(rozwiazanie(s+ds,lamstar)[1]-rozwiazanie(s-ds,lamstar)[1])/(2*ds)
    
return ret
end



@time lampos2(1,10)
  0.009549 seconds (144.32 k allocations: 3.936 MiB)
Out[10]:
-0.09663934541492836
In [11]:
function Fslam(s,lamstar)
return rozwiazanie(s,lamstar)[3]
end


function FslamS(lamstar)
return -1/(1-szift(lamstar,xpoczU))
end






function intelamstarU(lamstar)
 
return FslamS(lamstar)*(-1)*gamma*sqrt(delta)
end

function oneoverintelamstarU(lamstar)
return 1/intelamstarU(lamstar)
end



function bisection(f, a, b, tol)
    
    if sign(f(a)) == sign(f(b))
        error("function has the same sign at given endpoints")
    end

    mid = (a + b)/2

    while abs(f(mid)) > tol

        sign(f(mid)) == sign(f(a)) ? a=mid : b=mid
            
        mid = (a + b)/2

    end

    return mid
    
end


lamstar_zero1 = bisection(oneoverintelamstarU, -100, 0, 1e-12)
lamstar_zero2 = bisection(oneoverintelamstarU, 0, 50, 1e-12)


function calklamstarU(lamstar1)

    if lamstar_zero1<lamstar1<lamstar_zero2

    calk=quadgk(lamstar->intelamstarU(lamstar),0,lamstar1)[1]
        
   elseif lamstar1<lamstar_zero1
    
    calk=(quadgk(lamstar->intelamstarU(lamstar),0,lamstar_zero1+1e-7)[1] +
        quadgk(lamstar->intelamstarU(lamstar),lamstar_zero1-1e-7,lamstar1)[1] )
    
        
    elseif lamstar_zero2<lamstar1
    
    calk=(quadgk(lamstar->intelamstarU(lamstar),0,lamstar_zero2-1e-9)[1] +
        quadgk(lamstar->intelamstarU(lamstar),lamstar_zero2+1e-9,lamstar1)[1]  )
        
    end   
return calk
end





function ufinalSlam(lamstar)
return 0+ calklamstarU(lamstar) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end


function ufinalin(l, x0)
return 0+ calklamstarU(lamstar2(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end



### sprawdzenie calkUslam ## JEST OK

#function dcalkUslam(s1,lamstar)
#ds1=1e-3
#return (calkUslam(s1+ds1,lamstar)-calkUslam(s1-ds1,lamstar))/(2*ds1)
#end

#star1=-10
#Plots.plot(range(-10,0,length=500), [s1->dcalkUslam(s1,star1), s1->inteUslam(s1,star1)] )
#ufinalin(1,10)
#star1=5
#Plots.plot(range(-10,0,length=500), s1->ufinalSlam(s1, star1) )



#Plots.plot(range(-120,120,length=500), lamstar-> intelamstarU(lamstar) )
Out[11]:
ufinalin (generic function with 1 method)
In [12]:
Plots.plot(range(-120,20,length=500), lamstar->ufinalSlam(lamstar))
Out[12]:
In [13]:
xpoczV=10

function rozwiazanieV(s1,lamstar)

if s1<0    
    
function calyproblem(du, u, p, s)
du[1] = -derxpox0(u[1], u[2])
du[2] = gamma*sqrt(delta)
du[3] = u[3]*derxpox0polam(u[1], u[2])
end

u0 = [lamstar; xpoczV; 1/(1+szift(lamstar,xpoczV))]

sspan=(0,-s1 ) 

prob3 = ODEProblem(calyproblem, u0, sspan, abstol = 1e-5, reltol = 1e-5)
sol3 = solve(prob3)

return sol3(-s1)
        
    elseif s1>=0        
        
function calyproblem2(du, u, p, s)
du[1] = derxpox0(u[1], u[2])
du[2] = -gamma*sqrt(delta)
du[3] = -u[3]*derxpox0polam(u[1], u[2])
end

u02 = [lamstar; xpoczV; 1/(1+szift(lamstar,xpoczV))]

sspan2=(0,s1 ) 

prob32 = ODEProblem(calyproblem2, u02, sspan2, abstol = 1e-5, reltol = 1e-5)
sol32 = solve(prob32)

return sol32(s1)        
        
end
end


plotlyjs()
p=Plots.plot()
for i1 in range(-40,140,length=30)

    
Plots.plot!( s-> rozwiazanieV(s,i1)[2], s->rozwiazanieV(s,i1)[1], range(-20,20,length=500) )
end

display(p)
In [14]:
function lambda1V(s,lamstar)
return rozwiazanieV(s,lamstar)[1]
end

function sodx0V(x0)
return (xpoczV-x0)/(gamma*sqrt(delta))
end


function lamstar1V(x0, lambda)
    
s=sodx0V(x0)    
tab1=range(-550,950,length=1000)
tab2=lambda1V.(s,tab1)

itp=linear_interpolation(tab2,tab1)

return itp(lambda)
end


function lamstar2V(x0, lambda)
    
s=sodx0V(x0)    

function test(l)    
return lambda1V(s,l)-lambda
end     

ret=find_zero(test,-1150)

return ret
end

#lamstar1(1,lambda1(-5,21))

x001=0.05
lamstar1V(x001,lambda1(sodx0(x001),100))

function FfinalV(lam, x0)
    
if x0<xpoczV   
    
ret = rozwiazanieV(sodx0V(x0),lamstar2V(x0,lam))[3]
return ret
        
        
    elseif x0>=xpoczV      
return 1/(1+szift(lam,x0)) 
            
end 
        
        
end
Out[14]:
FfinalV (generic function with 1 method)
In [15]:
function derFpox0V(lam,x0)

function test(x0)
return FfinalV(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end


function derFpox02V(lam,x0)

function test(x0)
return FfinalV(lam, x0)
end
    
dx0=1e-6
    
return (test(x0+dx0)-test(x0))/dx0
end



function derF_PRAWOpolamV(lam, x0)
    
    
function test(lam)
return FfinalV(lam,x0)*derxpox0(lam, x0)
end
    
dlam=1e-6
    
return (test(lam+dlam)-test(lam))/dlam
end
Out[15]:
derF_PRAWOpolamV (generic function with 1 method)
In [16]:
function FslamV(s,lamstar)
return rozwiazanieV(s,lamstar)[3]
end


function FslamSV(lamstar)
return 1/(1+szift(lamstar,xpoczV))
end



function intelamstarV(lamstar)
s=0    
return FslamSV(lamstar)*gamma*sqrt(delta)
end


function oneoverintelamstarV(lamstar)
return 1/intelamstarV(lamstar)
end

#Plots.plot(range(-20,30,length=500), lamstar -> oneoverintelamstarV(lamstar))


Vlamstar_zero1 = bisection(oneoverintelamstarV, 0, 7, 1e-12)
Vlamstar_zero2 = bisection(oneoverintelamstarV, 6, 105, 1e-12)




function calklamstarV(lamstar1)

    if lamstar1<Vlamstar_zero1

    calk=quadgk(lamstar->intelamstarV(lamstar),0,lamstar1)[1]
        
    elseif Vlamstar_zero1<lamstar1<Vlamstar_zero2
    
    calk=(quadgk(lamstar->intelamstarV(lamstar),0,Vlamstar_zero1-1e-9)[1] +
        quadgk(lamstar->intelamstarV(lamstar),Vlamstar_zero1+1e-9,lamstar1)[1] )
    
        
    elseif Vlamstar_zero2<lamstar1
    
    calk=(quadgk(lamstar->intelamstarV(lamstar),0,Vlamstar_zero1-1e-9)[1] +
        quadgk(lamstar->intelamstarV(lamstar),Vlamstar_zero1+1e-9,Vlamstar_zero2-1e-7)[1] 
        +quadgk(lamstar->intelamstarV(lamstar),Vlamstar_zero2+1e-7,lamstar1)[1]  )
        
   end   
return calk
end

function vfinalSlam(lamstar)
return 0+ calklamstarV(lamstar) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end

function vfinalin(l, x0)
return 0+ calklamstarV(lamstar2V(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end



@time Plots.plot(range(-10,100,length=200), l->vfinalin(l, 10) )
102.284598 seconds (2.64 G allocations: 40.148 GiB, 3.79% gc time, 4.48% compilation time)
Out[16]:
In [17]:
##########################################
#####################################
######################################
######################################



function invUlam(lam,x0)
return 1/(tprim(lam)- 1/(1-szift(lam,x0))*xprim(lam, x0))
end


Plots.plot(range(-500,120,length=5000), lam->invUlam(lam,100) )

function invUlamzeros(x0)
    
    
function invUlamtest(lam)
return 1/(tprim(lam)- 1/(1-szift(lam,x0))*xprim(lam, x0))
end
    
return find_zeros(invUlamtest, -250,250)
end


function inteUlam(lam,x0)
return (tprim(lam)- 1/(1-szift(lam,x0))*xprim(lam, x0))
end

function calkUlam(lam,x0)
    

    
if isempty(invUlamzeros(x0))
 
calk = quadgk( l -> inteUlam(l, x0), 0, lam)[1]

    else  

        
 
        
        
        
        
        
    zero1=invUlamzeros(x0)[1]
zero2=invUlamzeros(x0)[2]
        
if zero1<0 && zero2>0        
        
        
        
        
        
if zero1<lam<zero2
calk = quadgk( l -> inteUlam(l, x0), 0, lam)[1]
        
elseif lam<zero1
#calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1+10e-4, rtol=10e-5 )[1]
#        + quadgk( l -> inteUlam(l, x0), zero1-10e-4, lam
#            , rtol=10e-5)[1] )

        
calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1+10e-4,  )[1]
        + quadgk( l -> inteUlam(l, x0), zero1-10e-4, lam
            )[1] )      
        
        
#elseif zero2<lam
#calk =( quadgk( l -> inteUlam(l, x0 ), 0, zero2-10e-6, rtol=10e-7)[1] 
#        + quadgk( l -> inteUlam(l, x0), zero2+10e-6, lam, 
#            rtol=10e-7)[1] )
 
        
elseif zero2<lam
calk =( quadgk( l -> inteUlam(l, x0 ), 0, zero2-10e-7)[1] 
        + quadgk( l -> inteUlam(l, x0), zero2+10e-7, lam, 
            )[1] )        
        
        
        
    end
            
        elseif zero1>0 && zero2>0            
          
            
            
            if lam<zero1
calk = quadgk( l -> inteUlam(l, x0), 0, lam)[1]
        
elseif zero1<lam<zero2
#calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1+10e-4, rtol=10e-5 )[1]
#        + quadgk( l -> inteUlam(l, x0), zero1-10e-4, lam
#            , rtol=10e-5)[1] )

        
calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1-10e-5,  )[1]
        + quadgk( l -> inteUlam(l, x0), zero1+10e-5, lam
            )[1] )      
        
        
#elseif zero2<lam
#calk =( quadgk( l -> inteUlam(l, x0 ), 0, zero2-10e-6, rtol=10e-7)[1] 
#        + quadgk( l -> inteUlam(l, x0), zero2+10e-6, lam, 
#            rtol=10e-7)[1] )
 
        
elseif zero2<lam
calk =(quadgk( l -> inteUlam(l, x0 ), 0, zero1-10e-5,  )[1]
        + quadgk( l -> inteUlam(l, x0), zero1+10e-5, zero2-10e-5)[1] +
                  quadgk( l -> inteUlam(l, x0), zero2+10e-5, lam)[1]  )   
            
            
            
            end          
            
            
        end           
    end

return calk
    

    
end

function xprimX0(x0)
lam1 = 0
function xxx1(x0)
lam1=0
beta0=Beta0(x0)
return (x0*(4+9*lam1^2 + 12*lam1*1/tan(beta0) +4*1/(tan(beta0))^2  )^(1/3))/(2^(2/3) * (1+ 1/(tan(beta0))^2 )^(1/3) )
end
  
dx0=10e-5
ret = (xxx1(x0+dx0)-xxx1(x0-dx0))/(2*dx0)
    
return ret
end

function xprimX000(x0)
return 1
end

function intenaX0(x0)
lam1 = 0
return -1/(1-szift(lam1, x0))*xprimX000(x0)
end

function oneoverintenaX0(x0)
return 1/intenaX0(x0)
end






function calkuodX0(x0)

if isempty(find_zeros(oneoverintenaX0, 0.1,1000))
        
calk = quadgk(x01 -> intenaX0(x01), 10, x0 )[1]

else        
    
zero1=find_zeros(oneoverintenaX0, 0.1,1000)[1]
zero2=find_zeros(oneoverintenaX0, 0.1,1000)[2]

if zero1<x0<zero2
calk = quadgk(x01 -> intenaX0(x01), 10, x0 )[1]
        
elseif x0<zero1
        
calk = quadgk(x01 -> intenaX0(x01), 10, zero1+10e-8 )[1] +  quadgk(x01 -> intenaX0(x01), zero1-10e-8, x0 )[1]

    elseif x0>zero2
calk=quadgk(x01 -> intenaX0(x01), 10, zero2-10e-6 )[1] +  quadgk(x01 -> intenaX0(x01), zero2+10e-6, x0 )[1]
    
end
    end        
        
return calk
end


function ufinal(lam, x0)
return calkUlam(lam,x0)+calkuodX0(x0)
end


@time ufinal(1,1)
  0.559610 seconds (721.59 k allocations: 31.450 MiB, 97.05% compilation time)
Out[17]:
-1.9976469928062972
In [18]:
function invVlam(lam,x0)
return 1/(tprim(lam)+ 1/(1+szift(lam,x0))*xprim(lam, x0))
end


#Plots.plot(range(-100,900,length=5000), lam->invVlam(lam,100) )



function invVlamzeros(x0)
    
    
function invVlamtest(lam)
return 1/(tprim(lam) + 1/(1+szift(lam,x0))*xprim(lam, x0))
end
    
return find_zeros(invVlamtest, -300,400)
end



function inteVlam(lam,x0)
return (tprim(lam)+ 1/(1+szift(lam,x0))*xprim(lam, x0))
end


function calkVlam(lam,x0)
    
if isempty(invVlamzeros(x0))
 
calk = quadgk( l -> inteVlam(l, x0), 0, lam)[1]

    else  
        
        
    zero1=invVlamzeros(x0)[1]
zero2=invVlamzeros(x0)[2]
    
if zero1<lam<zero2
calk = (quadgk( l -> inteVlam(l, x0 ), 0, zero1-10e-5)[1]
        + quadgk( l -> inteVlam(l, x0), zero1+10e-5, lam
            , rtol=10e-6)[1] )
        
elseif lam<zero1
        
calk = quadgk( l -> inteVlam(l, x0), 0, lam, rtol=10e-5 )[1]    
        
elseif zero2<lam
calk =(quadgk( l -> inteVlam(l, x0 ), 0, zero1-10e-5)[1]
        + quadgk( l -> inteVlam(l, x0), zero1+10e-5, zero2-10e-3)[1] +  quadgk( l -> inteVlam(l, x0), zero2+10e-3, lam)[1]  )
    
    end

return calk
    

    
end

end





function VintenaX0(x0)
lam1 = 0
return 1/(1+szift(lam1, x0))*xprimX000(x0)
end

function VoneoverintenaX0(x0)
return 1/VintenaX0(x0)
end


#Plots.plot(range(-1,30, length=1000), x0->VintenaX0(x0))


function VcalkuodX0(x0)

calk = quadgk(x0->VintenaX0(x0), 10, x0)[1]
return calk
end



function vfinal(lam, x0)
return calkVlam(lam,x0)+VcalkuodX0(x0)
end
Out[18]:
vfinal (generic function with 1 method)
In [19]:
######## PRE COMPUTING U I V 

#lamstarTAB=range(-100,100,length=5000)
#ufinalSlamTAB=ufinalSlam.(lamstarTAB)



#ufinalSlamINT=linear_interpolation(lamstarTAB,ufinalSlamTAB)


#function ufinalinINT(l, x0)
#return ufinalSlamINT(lamstar1(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
#end


####### V

#lamstarTAB=range(-100,100,length=5000)
#vfinalSlamTAB=vfinalSlam.(lamstarTAB)



#vfinalSlamINT=linear_interpolation(lamstarTAB,vfinalSlamTAB)



#function vfinalinINT(l, x0)
#return vfinalSlamINT(lamstar1(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
#end
In [20]:
kap=1/5


function atanufinal(lam, x0)
return 1/pi*atan(kap*ufinalin(lam,x0))
end


    function atanufinal2(lam, x0)
    return -1/pi*atan(kap*ufinalin(lam,x0))+1
    end

          function atanufinal3(lam, x0)
            return 1/pi*atan(kap*ufinalin(lam,x0))+2
            end


function atanvfinal(lam, x0)
return 1/pi*atan(kap*vfinalin(lam,x0))
end

    function atanvfinal2(lam, x0)
    return -1/pi*atan(kap*vfinalin(lam,x0))+1
    end

        function atanvfinal3(lam, x0)
        return 1/pi*atan(kap*vfinalin(lam,x0))+2
        end








################# TRAJEKTORIA X0=10


Ax0=10
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)

vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)



tab1_A = range(-200,ulamsin1-10e-5,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=50)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-7,length=50)
tab4_A = range(vlamsin1+10e-7, vlamsin2-10e-3,length=50)
tab5_A = range(vlamsin2+10e-3, 300,length=50)


tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamAa=tablamA


utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)

tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])

vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)

tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])



utrajA=linear_interpolation(tablamA, tabutabA );
vtrajA=linear_interpolation(tablamA, tabvtabA );



####  TRAJEKTORIA X0=6


Ax0=5
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)

vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)



tab1_A = range(-200,ulamsin1-10e-1,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=150)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-7,length=50)
tab4_A = range(vlamsin1+10e-7, vlamsin2-10e-3,length=50)
tab5_A = range(vlamsin2+10e-3, 300,length=50)


tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamBb=tablamA


utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)

tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])

vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)

tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])



utrajB=linear_interpolation(tablamA, tabutabA );
vtrajB=linear_interpolation(tablamA, tabvtabA );



################## TRAJEKTORIA X0=0.001


Ax0=0.001
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)

vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)



tab1_A = range(-400,ulamsin1-10e-3,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=150)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-7,length=50)
tab4_A = range(vlamsin1+10e-7, vlamsin2-10e-3,length=50)
tab5_A = range(vlamsin2+10e-3, 300,length=50)


tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamCc=tablamA


utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)

tabutabA =  vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])

vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)

tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])



utrajC=linear_interpolation(tablamA, tabutabA );
vtrajC=linear_interpolation(tablamA, tabvtabA );
In [21]:
############# TRAJEKTORIA X0=11


## UZYTE BYLE INNE U V 

function OUTatanufinal(lam, x0)
return 1/pi*atan(kap*ufinal(lam,x0))
end


    function OUTatanufinal2(lam, x0)
    return -1/pi*atan(kap*ufinal(lam,x0))+1
    end

          function OUTatanufinal3(lam, x0)
            return 1/pi*atan(kap*ufinal(lam,x0))+2
            end


function OUTatanvfinal(lam, x0)
return 1/pi*atan(kap*vfinal(lam,x0))
end

    function OUTatanvfinal2(lam, x0)
    return -1/pi*atan(kap*vfinal(lam,x0))+1
    end

        function OUTatanvfinal3(lam, x0)
        return 1/pi*atan(kap*vfinal(lam,x0))+2
        end




Ax0=11

Asin=invUlamzeros(Ax0)[1]
Asinb=invUlamzeros(Ax0)[2]
Asinc=invVlamzeros(Ax0)[1]
Asind=invVlamzeros(Ax0)[2]


tab1_A=range(-300,Asin-10e-1, length=250)
tab2_A=range(Asin+10e-1,Asinb-10e-4 , length=150)
tab3_A=range(Asinb+10e-4,Asinc-10e-4 , length=50)
tab4_A=range(Asinc+10e-4,Asind-10e-4 , length=50)
tab5_A=range(Asind+10e-4,300 , length=50)

tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamDd=tablamA

utab1_A = OUTatanufinal.(tab1_A, Ax0)
utab2_A = OUTatanufinal2.(tab2_A, Ax0)
utab3_A = OUTatanufinal3.(tab3_A, Ax0)
utab4_A = OUTatanufinal3.(tab4_A, Ax0)
utab5_A = OUTatanufinal3.(tab5_A, Ax0)

tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])



vtab1_A = OUTatanvfinal.(tab1_A, Ax0)
vtab2_A = OUTatanvfinal.(tab2_A, Ax0)
vtab3_A = OUTatanvfinal.(tab3_A, Ax0)
vtab4_A = OUTatanvfinal2.(tab4_A, Ax0)
vtab5_A = OUTatanvfinal3.(tab5_A, Ax0)

tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])


utrajD=linear_interpolation(tablamA, tabutabA );
vtrajD=linear_interpolation(tablamA, tabvtabA );



################ TRAJEKTORIA X0=15

Ax0=15

Asin=invUlamzeros(Ax0)[1]
Asinb=invUlamzeros(Ax0)[2]
Asinc=invVlamzeros(Ax0)[1]
Asind=invVlamzeros(Ax0)[2]


tab1_A=range(-300,Asin-10e-1, length=250)
tab2_A=range(Asin+10e-1,Asinb-10e-4 , length=150)
tab3_A=range(Asinb+10e-4,Asinc-10e-4 , length=50)
tab4_A=range(Asinc+10e-4,Asind-10e-4 , length=50)
tab5_A=range(Asind+10e-4,300 , length=50)

tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamEe=tablamA

utab1_A = OUTatanufinal.(tab1_A, Ax0)
utab2_A = OUTatanufinal2.(tab2_A, Ax0)
utab3_A = OUTatanufinal3.(tab3_A, Ax0)
utab4_A = OUTatanufinal3.(tab4_A, Ax0)
utab5_A = OUTatanufinal3.(tab5_A, Ax0)

tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])



vtab1_A = OUTatanvfinal.(tab1_A, Ax0)
vtab2_A = OUTatanvfinal.(tab2_A, Ax0)
vtab3_A = OUTatanvfinal.(tab3_A, Ax0)
vtab4_A = OUTatanvfinal2.(tab4_A, Ax0)
vtab5_A = OUTatanvfinal3.(tab5_A, Ax0)

tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])


utrajE=linear_interpolation(tablamA, tabutabA );
vtrajE=linear_interpolation(tablamA, tabvtabA );
In [22]:
####  TRAJEKTORIA X0=5


Ax0=5
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)

vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)



tab1_A = range(-300,ulamsin1-10e-1,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=150)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-10,length=150)
tab4_A = vcat([vlamsin1+10e-7] ,range(vlamsin1+8, vlamsin2-10e-3,length=50))
tab5_A = range(vlamsin2+10e-3, 300,length=50)


tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamBb=tablamA


utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)

tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])

vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)

tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])



utrajB=linear_interpolation(tablamA, tabutabA );
vtrajB=linear_interpolation(tablamA, tabvtabA );
In [23]:
[ulamsin1, ulamsin2, vlamsin1, vlamsin2]
Out[23]:
4-element Vector{Float64}:
 -110.63904313261818
   -5.236859665981674
    6.757238501613085
  116.048470561546
In [24]:
[ulamsin1, ulamsin2, vlamsin1, vlamsin2]
Out[24]:
4-element Vector{Float64}:
 -110.63904313261818
   -5.236859665981674
    6.757238501613085
  116.048470561546
In [25]:
function lamsin1(x0)
ulamsin1=lambda1(sodx0(x0) ,lamstar_zero1)
return ulamsin1
end

function lamsin2(x0)
ulamsin2=lambda1(sodx0(x0) ,lamstar_zero2)
return ulamsin2
end

function lamsin3(x0)
vlamsin1=lambda1V(sodx0V(x0) ,Vlamstar_zero1)
return vlamsin1
end

function lamsin4(x0)
vlamsin2=lambda1V(sodx0V(x0) ,Vlamstar_zero2)
return vlamsin2
end
Out[25]:
lamsin4 (generic function with 1 method)
In [26]:
#pv
function minusNzeros(x0)
function test(lam)
return 1-szift(lam, x0)
    end
    
    
lamtab1=range(-1000,1750,length=40000)
appalamtab=[]
    for i1 in lamtab1
    push!(appalamtab, test(i1))
end

miejsce=lamtab1[findmin(appalamtab)[2]]
    
    
return find_zeros(test, (miejsce-150,miejsce+150))
end





function plusNzeros(x0)
function test(lam)
return 1+szift(lam, x0)
    end
    
    
lamtab1=range(-1000,1750,length=40000)
appalamtab=[]
    for i1 in lamtab1
    push!(appalamtab, test(i1))
end

miejsce=lamtab1[findmin(appalamtab)[2]]
    
    
return find_zeros(test, (miejsce-150,miejsce+450))
end



st=2.97
Plots.plot(range(0.1,25,length=500), x0->lamsin1(x0), linewidth=2)
Plots.plot!(range(0.1,25,length=500), x0->lamsin2(x0), linewidth=2)
Plots.plot!(range(0.1,25,length=500), x0->lamsin3(x0), linewidth=2)
Plots.plot!(range(0.1,25,length=500), x0->lamsin4(x0), linewidth=2)


Plots.plot!(range(st,25,length=500), x0->minusNzeros(x0)[1], color=:black)
Plots.plot!(range(st,25,length=500), x0->minusNzeros(x0)[2], color=:black)
Plots.plot!(range(st,25,length=500), x0->plusNzeros(x0)[1], color=:black)
Plots.plot!(range(st,25,length=500), x0->plusNzeros(x0)[2], color=:black, xlims=(0,20))
Out[26]:
In [27]:
function tescik(x0)
return minusNzeros(x0)[1] - lamsin2(x0)
end

#tescik(5.25)

przec1 = find_zero(tescik, (st,10) )
Out[27]:
5.061379984358704
In [28]:
#### APPARENT HORIZONS

















function app_v_dol1(x0)
return atanvfinal(minusNzeros(x0)[1],x0)
end

function app_v_dol2(x0)
return atanvfinal(minusNzeros(x0)[2],x0)
end


function app_u_dol1(x0)
return atanufinal2(minusNzeros(x0)[1],x0)
end

function app_u_dol2(x0)
return atanufinal2(minusNzeros(x0)[2],x0)
end



function app_u_dol1W(x0)
return atanufinal3(minusNzeros(x0)[1],x0)
end

function app_u_dol2W(x0)
return atanufinal3(minusNzeros(x0)[2],x0)
end






### APP DOL INSIDE
x0ran=range(przec1,9.999,length=150)
app_vtab1=app_v_dol1.(x0ran)
app_utab1=app_u_dol1.(x0ran)

linapp_v = linear_interpolation(x0ran, app_vtab1)
linapp_u = linear_interpolation(x0ran, app_utab1)


Plots.plot( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=:magenta )



x0ran2 = range(st, przec1, length=150)

app_vtab1_przec = app_v_dol1.(x0ran2)
app_utab1_przec = app_u_dol1W.(x0ran2)

linapp_v_przec=linear_interpolation(x0ran2,app_vtab1_przec)
linapp_u_przec=linear_interpolation(x0ran2, app_utab1_przec)




Plots.plot!( x0->linapp_v_przec(x0) - linapp_u_przec(x0), x0->linapp_v_przec(x0) + linapp_u_przec(x0), x0ran2, color=:magenta )

Plots.plot!( [linapp_v(przec1) - linapp_u(przec1),linapp_v_przec(przec1) - linapp_u_przec(przec1) ],
             [linapp_v(przec1) + linapp_u(przec1),linapp_v_przec(przec1) + linapp_u_przec(przec1) ], color=:magenta )


x0ran3= range(st, 10, length=150)

app_vtab2=app_v_dol2.(x0ran3)
app_utab2=app_u_dol2W.(x0ran3)


linapp_v2 = linear_interpolation(x0ran3, app_vtab2)
linapp_u2 = linear_interpolation(x0ran3, app_utab2)


Plots.plot!( x0->linapp_v2(x0) - linapp_u2(x0), x0->linapp_v2(x0) + linapp_u2(x0), x0ran3, color=:magenta )


Plots.plot!( [linapp_v_przec(st) - linapp_u_przec(st), linapp_v2(st) - linapp_u2(st)],
[linapp_v_przec(st) + linapp_u_przec(st), linapp_v2(st) + linapp_u2(st)], color=:magenta )



Plots.plot!(  [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , 0.5 -0.5  ], [OUTatanvfinal(lamsin1(10) ,10) + 0.5 , 0.5 +0.5  ], color=:magenta   )

### APP DOL OUTSIDE

#x0ran=range(11,50,length=150)


#app_vtab1=app_v_dol1.(x0ran)
#app_utab1=app_u_dol1.(x0ran)

#linapp_vOUT = linear_interpolation(x0ran, app_vtab1)
#linapp_uOUT = linear_interpolation(x0ran, app_utab1)

#Plots.plot!( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=:blue )
Out[28]:
In [29]:
### APP GORNY

function tescik2(x0)
return plusNzeros(x0)[2] - lamsin3(x0)
end

#tescik(5.25)

przec2 = find_zero(tescik2, (st,10) )


function upAPPv(x0)
return atanvfinal(plusNzeros(x0)[1],x0)
end

function uppAPPu(x0)
return atanufinal3(plusNzeros(x0)[1],x0)
end

function upAPPvB(x0)
return atanvfinal(plusNzeros(x0)[2],x0)
end

function uppAPPuB(x0)
return atanufinal3(plusNzeros(x0)[2],x0)
end


function upAPPvBW(x0)
return atanvfinal2(plusNzeros(x0)[2],x0)
end

function uppAPPuBW(x0)
return atanufinal3(plusNzeros(x0)[2],x0)
end





up_x0ran=range(st,9.999, length=150)

uptabav1= upAPPv.(up_x0ran)
uptabau1= uppAPPu.(up_x0ran)


linupav1 = linear_interpolation(up_x0ran, uptabav1)
linupau1 = linear_interpolation(up_x0ran, uptabau1)


## APP wyzej


up_x0ran2=range(st, przec2, length=150)

uptabavB= upAPPvB.(up_x0ran2)
uptabauB= uppAPPuB.(up_x0ran2)


linupavB = linear_interpolation(up_x0ran2, uptabavB)
linupauB = linear_interpolation(up_x0ran2, uptabauB)



## APP PO PRZECIECIU 

up_x0ran3=vcat(range(przec2,przec2+1e-3 , length=75), range(4.1,9.999,length=75))

uptabavBW= upAPPvBW.(up_x0ran3)
uptabauBW= uppAPPuBW.(up_x0ran3)


linupavBW = linear_interpolation(up_x0ran3, uptabavBW)
linupauBW = linear_interpolation(up_x0ran3, uptabauBW)




Plots.plot( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=:darkviolet )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=:darkviolet )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=:darkviolet )


Plots.plot!( [linupav1(st)-linupau1(st) ,linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ], color=:darkviolet  )


Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
    color=:darkviolet  )
Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
    [linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
    color=:darkviolet  )
Out[29]:
In [30]:
Plots.plot( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=:darkviolet )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=:darkviolet )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=:darkviolet )


Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
    color=:darkviolet  )

Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
    [linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
    color=:darkviolet  )
Out[30]:
In [31]:
## OBA APPARENT

Plots.plot( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=:magenta )
Plots.plot!( x0->linapp_v_przec(x0) - linapp_u_przec(x0), x0->linapp_v_przec(x0) + linapp_u_przec(x0), x0ran2, color=:magenta )
Plots.plot!( x0->linapp_v2(x0) - linapp_u2(x0), x0->linapp_v2(x0) + linapp_u2(x0), x0ran3, color=:magenta )



Plots.plot!( [linapp_v(przec1) - linapp_u(przec1),linapp_v_przec(przec1) - linapp_u_przec(przec1) ],
             [linapp_v(przec1) + linapp_u(przec1),linapp_v_przec(przec1) + linapp_u_przec(przec1) ], color=:magenta )


Plots.plot!( [linapp_v_przec(st) - linapp_u_przec(st), linapp_v2(st) - linapp_u2(st)],
[linapp_v_przec(st) + linapp_u_przec(st), linapp_v2(st) + linapp_u2(st)], color=:magenta )

Plots.plot!(  [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , 0.5 -0.5  ], [OUTatanvfinal(lamsin1(10) ,10) + 0.5 , 0.5 +0.5  ], color=:magenta   )
Plots.plot!(  [OUTatanvfinal(lamsin2(10) ,10) - 1.5 , 0.5 -1.5  ], [OUTatanvfinal(lamsin2(10) ,10) + 1.5 , 0.5 +1.5  ], color=:magenta   )

Plots.plot!([OUTatanvfinal(lamsin2(10) ,10) - 1.5 , linapp_v2(9.999) - linapp_u2(9.999) ], 
    [OUTatanvfinal(lamsin2(10) ,10) + 1.5 , linapp_v2(9.999) + linapp_u2(9.999) ], color=:magenta)



Plots.plot!(  [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , linapp_v(9.999) - linapp_u(9.999) ], 
    [OUTatanvfinal(lamsin1(10) ,10) + 0.5 ,  linapp_v(9.999) + linapp_u(9.999)  ], color=:magenta   )



## DOLNY
Plots.plot!( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=:darkviolet )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=:darkviolet )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=:darkviolet )


Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
    color=:darkviolet  )

Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
    [linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
    color=:darkviolet  )
Out[31]:
In [32]:
## TIMELIKE SINGULARITY

function trajx(T)
return x(T/(gamma*sqrt(delta)),10)
end

M= 4/3*pi*xb^3*rho0
ff = x -> 1 - 2*G*M/x + 4*G^2*M^2*gamma^2*delta/x^4


xs1=find_zeros(ff,(0,150) )[1]
xs2=find_zeros(ff,(0,150) )[2]


inte =x-> 1/(ff(x) )


xszift = x-> sqrt(2*G*M*(x^3-2*G*M*gamma^2*delta)/x^4)


function calkX(x)

if x<xs1
ret = quadgk( x1->inte(x1),  1e-15, x)[1]

elseif xs1<x<xs2
ret = quadgk( x1->inte(x1),  0, xs1-1e-7)[1] + quadgk( x1->inte(x1),  xs1+1e-7, x)[1]

elseif xs2<x

ret = quadgk( x1->inte(x1), 0, xs1-1e-7)[1] + quadgk( x1->inte(x1), xs1+1e-12, xs2-1e-7)[1] + quadgk( x1->inte(x1), xs2+1e-7, x)[1] 
end
    
return ret
end

intet = T-> 1/ff(trajx(T))

tsin1 = find_zeros(intet, (-100,100) )[1]
tsin2 = find_zeros(intet, (-100,100) )[2]
tsin3 = find_zeros(intet, (-100,100) )[3]
tsin4 = find_zeros(intet, (-100,100) )[4]




function calkt(T)

if T<tsin1
ret= quadgk( T1->intet(T1), 0, tsin1+1e-7 )[1] + quadgk( T1->intet(T1), tsin1-1e-7, T )[1]




elseif tsin1<T<tsin2
ret =  quadgk( T1->intet(T1), 0, T )[1]

elseif tsin2<T<tsin3
ret= quadgk( T1->intet(T1), 0, tsin2-1e-7 )[1] + quadgk( T1->intet(T1), tsin2+1e-7, T)[1]

elseif tsin3<T<tsin4
ret = (quadgk( T1->intet(T1), 0, tsin2-1e-7 )[1] + quadgk( T1->intet(T1), tsin2+1e-7, tsin3-1e-7 )[1]
            +quadgk( T1->intet(T1), tsin3+1e-7, T )[1] )

elseif tsin4<T
ret = (quadgk( T1->intet(T1), 0, tsin2-1e-7 )[1] + quadgk( T1->intet(T1), tsin2+1e-7, tsin3-1e-7 )[1]+
            quadgk( T1->intet(T1), tsin3+1e-7, tsin4-1e-7 )[1]
    +quadgk( T1->intet(T1), tsin4+1e-7, T)[1] )
end

return ret
end


function usin(t,x)
return -1/pi*atan( calkt(t) - calkX(x)  ) +1
end


function vsin(t,x)
return -1/pi*atan( calkt(t) + calkX(x)  ) +1
end

x1=0.5
Plots.plot( t-> vsin(t,x1)-usin(t,x1), t-> vsin(t,x1)+usin(t,x1), range(-10,10,length=100), color=:black )
Out[32]:
In [33]:
######### POWIERZCHNIA POCZATKOWA

in_przec1=find_zeros(lamsin2,(0,10))[1]
in_przec2=61.5558

#Plots.plot(x0->atanvfinal(0,x0)-atanufinal3(0,x0), x0->atanvfinal(0,x0)+atanufinal3(0,x0), range(1e-5,in_przec1,length=100), color=:black)
#Plots.plot!(x0->atanvfinal(0,x0)-atanufinal2(0,x0), x0->atanvfinal(0,x0)+atanufinal2(0,x0), range(in_przec1,10,length=100), color=:black)

#Plots.plot!(x0->OUTatanvfinal(0,x0)-OUTatanufinal2(0,x0), x0->OUTatanvfinal(0,x0)+OUTatanufinal2(0,x0), range(10,in_przec2-1e-1,length=100),
#   color=:black)
#Plots.plot!(x0->OUTatanvfinal(0,x0)-OUTatanufinal(0,x0), x0->OUTatanvfinal(0,x0)+OUTatanufinal(0,x0),
#range(in_przec2+1e-1,120,length=100), color=:black)

intabx01=range(1e-5,in_przec1,length=100)
intabx02=range(in_przec1,10-1e-15,length=100)
intabx03=range(10+1e-15,in_przec2-1e-1,length=100)
intabx04=range(in_przec2+1e-1,120,length=100)

intabx0=vcat(intabx01,intabx02,intabx03,intabx04, [1000])


intabv1=atanvfinal.(0,intabx01)
intabv2=atanvfinal.(0,intabx02)
intabv3=OUTatanvfinal.(0,intabx03)
intabv4=OUTatanvfinal.(0,intabx04)


intabu1=atanufinal3.(0,intabx01)
intabu2=atanufinal2.(0,intabx02)
intabu3=OUTatanufinal2.(0,intabx03)
intabu4=OUTatanufinal.(0,intabx04)

intabv=vcat(intabv1,intabv2,intabv3,intabv4, [0.5])
intabu=vcat(intabu1,intabu2,intabu3,intabu4, [-0.5])

innv=linear_interpolation( intabx0, intabv)
innu=linear_interpolation( intabx0, intabu)

Plots.plot( x0-> innv(x0)-innu(x0), x0-> innv(x0)+innu(x0), intabx0, color=:black)
┌ Warning: Duplicated knots were deduplicated. Use Interpolations.deduplicate_knots!(knots) explicitly to avoid this warning.
│   k1 =
│    401-element Vector{Float64}:
│        1.0e-5
│        0.07401944551836069
│        0.14802889103672137
│        0.22203833655508207
│        0.29604778207344273
│        ⋮
│      118.23199393939394
│      118.82132929292929
│      119.41066464646465
│      120.0
│     1000.0
└ @ Interpolations C:\Users\mihalbobula\.julia\packages\Interpolations\91PhN\src\gridded\gridded.jl:77
Out[33]:
In [37]:
############ DIAGRAM

#plotlyjs()
gr()



######## APPPARENT DOLNY
lwa=1.4
cap=:magenta

Plots.plot( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=cap,
    label=L"\mathrm{first} \, \mathrm{AH}", lw=lwa )
Plots.plot!( x0->linapp_v_przec(x0) - linapp_u_przec(x0), x0->linapp_v_przec(x0) + linapp_u_przec(x0), x0ran2, color=cap, label=false, lw=lwa )
Plots.plot!( x0->linapp_v2(x0) - linapp_u2(x0), x0->linapp_v2(x0) + linapp_u2(x0), x0ran3, color=cap, label=false, lw=lwa )



Plots.plot!( [linapp_v(przec1) - linapp_u(przec1),linapp_v_przec(przec1) - linapp_u_przec(przec1) ],
             [linapp_v(przec1) + linapp_u(przec1),linapp_v_przec(przec1) + linapp_u_przec(przec1) ], color=cap, label=false, lw=lwa )


Plots.plot!( [linapp_v_przec(st) - linapp_u_przec(st), linapp_v2(st) - linapp_u2(st)],
[linapp_v_przec(st) + linapp_u_przec(st), linapp_v2(st) + linapp_u2(st)], color=cap, label=false, lw=lwa )

Plots.plot!(  [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , 0.5 -0.5  ], [OUTatanvfinal(lamsin1(10) ,10) + 0.5 , 0.5 +0.5  ],
    color=cap, label=false , lw=lwa  )
Plots.plot!(  [OUTatanvfinal(lamsin2(10) ,10) - 1.5 , 0.5 -1.5  ], [OUTatanvfinal(lamsin2(10) ,10) + 1.5 , 0.5 +1.5  ], color=cap, 
    label=false, lw=lwa   )

Plots.plot!([OUTatanvfinal(lamsin2(10) ,10) - 1.5 , linapp_v2(9.999) - linapp_u2(9.999) ], 
    [OUTatanvfinal(lamsin2(10) ,10) + 1.5 , linapp_v2(9.999) + linapp_u2(9.999) ], color=cap, label=false,  lw=lwa)



Plots.plot!(  [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , linapp_v(9.999) - linapp_u(9.999) ], 
    [OUTatanvfinal(lamsin1(10) ,10) + 0.5 ,  linapp_v(9.999) + linapp_u(9.999)  ], color=cap, label=false , lw=lwa  )



#### APPARENT GORNY

cap2=:cyan3

Plots.plot!( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=cap2,
    label=L"\mathrm{second} \, \mathrm{AH}", lw=lwa  )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=cap2, label=false, lw=lwa )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=cap2, label=false, lw=lwa )


Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
    color=cap2, label=false, lw=lwa  )

Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
    [linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
    color=cap2, label=false, lw=lwa  )


Plots.plot!( [0.5-OUTatanufinal3(lamsin3(10),10)  , 0.5-1.5  ],[0.5+OUTatanufinal3(lamsin3(10),10)  , 0.5+1.5  ],
    color=cap2, label=false, lw=lwa )
Plots.plot!( [1.5-OUTatanufinal3(lamsin4(10),10)  , 1.5-1.5  ],[1.5+OUTatanufinal3(lamsin4(10),10)  , 1.5+1.5  ],color=cap2,
    label=false, lw=lwa  )

Plots.plot!( [0.5-OUTatanufinal3(lamsin3(10),10)  , linupav1(9.999)-linupau1(9.999) ],
              [0.5+OUTatanufinal3(lamsin3(10),10)  , linupav1(9.999)+linupau1(9.999) ], color=cap2, label=false, lw=lwa  )


Plots.plot!( [1.5-OUTatanufinal3(lamsin4(10),10)  , linupavBW(9.999)-linupauBW(9.999)  ],
    [1.5+OUTatanufinal3(lamsin4(10),10)  , linupavBW(9.999)+linupauBW(9.999)  ],color=cap2, label=false, lw=lwa  )


### HORYZONTY ZEWNETRZNE

Plots.plot!( [0.5-0.5, 0.5-1.5 ], [0.5+0.5, 0.5+1.5 ], color=:grey, label=L"x_- \, \mathrm{horizons}"  )
Plots.plot!(  [0.5-1.5, 1.5-1.5 ] , [0.5+1.5, 1.5+1.5 ], color=:grey, label=false )

######### NULL INFINITIES
nl=:darkgreen

Plots.plot!( [-0.5+0.5, 0.5+0.5], [-0.5-0.5, 0.5-0.5], color=nl, label=L"\mathrm{null} \, \mathrm{infinities}"  )
Plots.plot!( [0.5-(-0.5) ,0.5-(0.5) ],  [0.5+(-0.5) ,0.5+(0.5) ], color=nl, label=false )

Plots.plot!( [1.5-1.5, 2.5-1.5], [1.5+1.5, 2.5+1.5], color=nl, label=false    )
Plots.plot!( [ 2.5-1.5, 2.5-2.5  ],  [ 2.5+1.5, 2.5+2.5  ], color=nl, label=false )


### TIMELIKE SINGULARITY
x1=0
Plots.plot!( t-> vsin(t,x1)-usin(t,x1), t-> vsin(t,x1)+usin(t,x1), range(-15000,15000,length=4), color=:red2, ls=:dot,linewidth=1,
label=L"{\mathrm{singularity}}")



######## POW POCZATKOWA

Plots.plot!( x0-> innv(x0)-innu(x0), x0-> innv(x0)+innu(x0), intabx0, color=:black,
    label=L"\mathrm{initial} \, \mathrm{surface}", ls=:dashdot  )




######## TRAJEKTORIE
lw1=0.6


Plots.plot!( l -> vtrajB(l) - utrajB(l),l -> vtrajB(l) + utrajB(l),tablamBb , aspectratio=1, color=:blue, label=L"x_0=5", lw=lw1  )
Plots.plot!( l -> vtrajC(l) - utrajC(l),l -> vtrajC(l) + utrajC(l),tablamCc , aspectratio=1, color=:green, label=L"x_0=0", lw=lw1  )
Plots.plot!( l -> vtrajD(l) - utrajD(l),l -> vtrajD(l) + utrajD(l),tablamDd , aspectratio=1, color=:orange, 
        linewidth=1, linestyle=:dash, label=L"x_0=11", lw=lw1)


Plots.plot!( l -> vtrajE(l) - utrajE(l),l -> vtrajE(l) + utrajE(l),tablamEe , aspectratio=1, color=:brown, 
        linewidth=1, linestyle=:dash, label=L"x_0=15", lw=lw1)


Plots.plot!( l -> vtrajA(l) - utrajA(l),l -> vtrajA(l) + utrajA(l), tablamAa, aspectratio=1, color=:red, label=L"x_0=10", lw=lw1  )
Out[37]:
In [38]:
###### ATRYBUTY

Plots.plot!(xlabel=L"\tilde{v} - \tilde{u}",
    ylabel=L"\tilde{v} + \tilde{u}")
Plots.plot!(xlims=(-4.25,3.75), ylims=(-1.25,5.25), grid=false, xticks=-3:4:1, yticks=-1:6:5, legendfontsize=6,
xtickfontsize=6,ytickfontsize=6, xguidefontsize=10,yguidefontsize=10 )



######### LENS

lens!([-2.2, -1.7], [1.6, 2.1];  inset=(1, bbox(0.1, 0.2,0.21, 0.21)), 
       subplot=2, ticks=false, framestyle=:box, ls=:dash, lw=0.5 , lc=:grey, aspectratio=1)



lens!([-1.1,-0.6 ], [-0.1, 0.4];  inset=(1, bbox(0.15, 0.72,0.21, 0.21)), 
       subplot=2, ticks=false, framestyle=:box, ls=:dash, lw=0.5 , lc=:grey, aspectratio=1)
Out[38]:
In [39]:
Plots.savefig( "dagINNY.svg")
Out[39]:
"C:\\Users\\mihalbobula\\Desktop\\ew_final_przec\\dagINNY.svg"
In [ ]: